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Abstract 

Red blood cells (RBCs) are an essential component of blood. A method to include the particulate nature of 
blood is introduced here with the goal of studying circulation in large-scale realistic vessels. The method uses a 
combination of the Lattice Boltzmann method (LBM) to account for the plasma motion, and a modified Molecular 
Dynamics scheme for the cellular motion. Numerical results illustrate the quality of the model in reproducing 

^-^nown rheological properties of blood as much as revealing the effect of RBC structuring on the wall shear stress, 

^3\fith consequences on the development of cardiovascular diseases. 
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I. INTRODUCTION 

The study of blood flows in physiological vessels, named hemodynamics, is an active fleld of research 
from a fundamental point of view and for the strategic medical implications related to cardiovascular 
diseases. Atherosclerosis, for instance, is the leading cause of death in western countries and its societal 
impact is constantly increasing [? ]. Atherogenesis is the formation of plaques within the endothelium, 
the tissue forming the arterial walls, and is triggered by low levels of shear stress at the wall region. The 
biomechanical origins of atherogenesis relate to the disturbed flow patterns encountered in close proximity 
of the endothelium [? ]. 

Hemodynamics entails several physical and biochemical levels and understanding blood flows in complex 
geometries, from the large-scale coronary arteries to microcapillaries, has attracted considerable attention 
over centuries, starting from the early investigations by Leonardo da Vinci [? ]. Blood exhibits both 
visco-elastic and thyxotropic response and, most importantly, shear-thinning, occurring at shear rates 
commonly encountered in flow conditions corresponding to large-scale arteries, such as in the coronary 
system. 

In the last decades, computer simulation has provided new understanding on the flow patterns of blood 
modeled as a continuum. In physiological condition, however, blood presents a large volume fraction of red 
blood cells (RBC), with hematocrit level H ranging between 35 and 50%. Recent computational models 
have been put forward to reproduce blood from a bottom-up perspective, that is, by including plasma 
and individual red blood cells, either as flexible bodies [???????] or as rigid ellipsoids [? ]. 
A detailed representation of RBCs is crucial to study microcirculation in capillaries, a situation where 
shape, deformability and near-fleld hydrodynamic response of single red blood cells need to be accurately 
accounted for. Conversely, not much effort has been devoted to modeling red blood cells in situations 
where the far-fleld hydrodynamics and the global rheological response of blood are taken into account, two 
aspects that can have strong impact on large-scale flows. Such possibility would enable to study blood flows 
in large cardiovascular networks, entailing from the red blood cell level up to heart-size geometries, with 
a computational effort that represents a trade-off between physical fldelity and computational feasibility. 

In this paper we address the issue of designing a robust physical representation of blood by a model 
that can be used in large-scale coronary arteries. Therefore, our model accounts for the role played by 
RBCs in the non-Newtonian hemorheology together with accounting for cell crowding in proximity of the 
vessel walls. This aspect is expected to have dramatic consequences on the distribution of endothelial 
shear stress, due to the inhomogeneous distribution of red blood cells in proximity of the arterial wall. 

Red blood cells are globules that present a biconcave discoidal form, and a soft membrane that encloses 
a high-viscosity liquid made of hemoglobin. In presence of an external shear fleld, red blood cells present 
both a solid-like tumbling and a vesicular motion with the attendant sliding of the membrane, the so-called 



tank treading motion. As a result, RBC encodes both rotational and orientational responses that deeply 
and non-trivially modulate blood rheology. We account for these contributions and neglect memory-related 
response of the cell membrane arising from the rearrangement of the cytoskeleton. In addition, we take 
into account some degree of cell softness for the RBC-RBC collisional properties and disregard global 
shape rearrangements. 

The proposed model focuses on three independent components: the far-field hydrodynamic interaction 
of a RBC in a plasma solvent, the raise of viscosity of the suspension with the hematocrit level and 
the many-body collisional contributions to viscosity. These three critical ingredients conspire to produce 
large-scale hemorheology and the local structuring of RBCs. 

The paper is organized as follows. In Section |TT] the computational model is described by focusing 
on blood plasma at first, and then incrementally introducing the RBC hydrodynamic interactions (roto- 
translational and tank treading coupling), the modulation of viscosity of the suspension and the RBC-RBC 



mechanical forces. Section III illustrates the numerical validation of the model and presents results on 



RBC-rich blood flows in a realistic vessel. 



II. COMPUTATIONAL MODEL 



A. Blood plasma 

More than 99 % in volume of blood is composed by plasma and red blood cells, where the latter can 
account up to 50 % in volume of the suspension. Plasma is a water-like Newtonian fluid that can be 
modeled as a continuum, with the Navier-Stokes equations governing the macro dynamics. 

For the purpose of reproducing the plasma evolution in high complex geometries (and vanishing hema- 
tocrit) together with the inclusion of red blood cells as suspended bodies, a convenient route is to model 
plasma at kinetic level via the Lattice Boltzmann (LB) method [? ]. The LB method deals with the 
evolution of the distribution function /(x, v, t), describing the probability to find a plasma particle at site 
X, velocity v and time t, and recovers Newtonian rheology in the limiting macroscopic space/time scale. 
By representing the distribution function in discretized form over a cartesian mesh x, /(x, v, t) — )■ /p(x, t), 
where the subscript p labels a set of discrete speeds Cp connecting mesh points to mesh neighbors, the 
evolution dynamics over a timestep h reads 

fp{^ + hcp,t + h) = /;(x,t) (1) 

with /* (x, t) being the post-coUisional population, 

/; = (1 - !^)fp + 1/;^ + hAff'^^ (2) 



T is a characteristic relaxation time and /^"^ the Maxwelhan equihbriuni expressed as a second-order 
low-Mach expansion in the fluid velocity u, 



fp'^ = WpP 



U • C„ (u • C„)^ — C^M^ 

1 + ^ + ^ ^^ 



(3) 



C2 2c4 

Eqs. ([ljl2j) encode the effect of streaming, that is, the kinematic motion of free particles along straight 
trajectories, together with the solvent-solvent and the solvent-solute "molecular" collisions. The solvent- 
solvent collisions are included via the so-called BGK relaxation kernel, expressed by the fact that the 
distribution strives to reach the local equilibrium over the timescale T. The plasma kinematic viscosity u 
relates to T via 

u = c\r-h/2) 

where c is the plasma sound speed. 

In this work, we employ the D3Q19 lattice scheme, where one has c = 1/v^ , and Wp stands 
for a set of normalized weights with p = 0,...,18, being equal to Wp = 1/3 for the population cor- 
responding to the null discrete speed Cq = (0,0,0), Wp = 1/18 for the ones connecting first mesh 
neighbors Ci e = (iljOjO)) (0, ±1,0), (0,0, ±1), and Wp = 1/36 for second mesh neighbors, Cy . , ig = 
(±1,±1,0), (±1,0,±1), (0,±1,±1). 

The term A/^'''^^ accounts for the presence of suspended RBC that act as body forces on the plasma 
in a hydrokinetic way, as detailed out in the next section. It has the following expression 



^fp' = hwpp 



G ■ Cp (G ■ Cp)(u ■ Cp) — c^G ■ u 



(4) 



c2 2c4 

where G is the local body-fluid coupling. Eq. (Il]) is a representation of the body force consistent with the 
second-order Hermite expansion of the distribution /(x, v,t) and BGK collisional kernel [? ]. 

Knowledge of the discrete populations fp allows to compute the local plasma density p, speed u and 
momentum- flux tensor P, by a direct summation over the populations 

P = Y.fv (5) 

p 

pu = ^ fpCp + -pG (6) 

p 

P = ^fpCpCp (7) 

p 
where eq. (Ig]) ensures second order space/time accuracy of the LB algorithm [? ]. 

The diagonal component of the momentum-fiux tensor gives the plasma pressure, while the off-diagonal 

terms give the shear stress cr, the latter being related to the non-equilibrium component of the populations, 

V 



up {du + dn') = ^Y1 ^P^P (fp - fp") 



On the other hand, due to the lack of a local kinetic definition of the antisymmetric component of the 
displacement tensor du and the fluid vorticity, these terms are evaluated via finite-differences. 



B. Red blood cells 



At first, red blood cells are introduced by modeling the hydrodynamic interaction of a red blood cell 
with the surrounding plasma solvent and next, by considering the collisional interaction of different red 
blood cells among themselves. We take an oblate ellipsoid as a good approximation to the hydrodynamic 
shape of a RBC 

A RBC is a body having mass M, position Rj, velocity Vj, angular velocity 17^, and instantaneous 
orientation given by the matrix 



Q^ 



/ 



^x,i '^x,i 9x,^ 
^y,i ^y,i 9y,i 



\ 



(9) 



z,i yz,i I 



where fij, tj. 



are orthogonal unit vectors, such that Qf Qj = 1. The orthogonal matrix Qj transforms 
between the body and the laboratory frame via v' = QjV, where the primed and unprimed symbols stand 
for laboratory and body frames, respectively. The tensor of inertia, Ij, is diagonal in the body frame and 
transforms to the laboratory frame according to l'^ = QjIjQ^ . In the sequel, we shall drop the prime 
symbol to ease the notation and implicitly mean that the translational motion is handled in the laboratory 
frame, where the rotational motion is handled in the body frame. 

We collectively denote the roto-translational state by the symbol Fj = (Rj, Qj, Vj, fij). Let us now 
introduce an auxiliary function to account for the shape and orientation of the suspended RBC. We choose 
the following expression [? ] 

5(x,Q,)= J] 54(Q,x)J (10) 



a=x,y,z 



with 



I (5 - 4:\yJU - \/l + Sbal/^a - 16l/2/e) 



^aiVa) = < 



I 3 - 4|t/J/e. - v/-7 + 24|t/,|/a-16|/2/e^ 







\ya/^a\ < 0.5 

0.5 < \yJU < 1 

\ya\/^a > 1 



'11^ 



and C,a being a set of three integers, one for each cartesian component a = x,y, z, representing the ellipsoidal 
radii in the three principal directions. The shape function has compact support and for C,x = C,y = ^z = '^ 
generates a spherically symmetric diffused particle with a support extending over 64 mesh points. Given 
that we handle one RBC via a single shape function, the computational cost of a suspended RBC is 
proportional to the size of the support in the three cartesian directions. 

The shape function has two important properties, it is normalized when summed over the cartesian 
mesh points x [? ] 

J]5(x-T) = l (12) 







f-JI^#^JIh 



Figure 1: The two types of motion considered in modeling a RBC in a linear shear field, the solid-like, tumbling 
or flipping coin motion (upper panel), and the vesicular, tank-treading motion (lower panel), where two material 
points on the RBC membrane move at fixed body orientation. 



for any continuous displacement T, and obeys the property 



5^(x„-T,)9^5(x-T) = -5, 



al3 



(13) 



The translational response of the suspended body is designed according to the RBC-fluid exchange 
kernel 

0(x, r,) = -7t5(x - R„ Q,,) [V, - u(x)] = -7t5, (V^ - u) (14) 

where Jt is a translational coupling coefficient and where the short-hand notation 6i = 5(x — Rj, Qj) has 
been introduced. 

The body rotational response has different origins and can be analyzed by considering the general 
decomposition of the deformation tensor in terms of purely elongational and rotational terms 

du = e + p 

where e = ^{du + du^) is the symmetric rate of strain tensor, related to the dissipative character of the 
flow, and p = ^{du — du^) is the antisymmetric vorticity tensor, which bears the conservative component 
of the flow and is related to the vorticity vector lo = d x u = e : p, where e is the Levi-Civita tensor [? ]. 
The rotational component of the deformation tensor gives rise to a solid-like tumbling motion, where the 
rotational and elongational one give rise to the vesicular, tank treading motion, as illustrated in Fig. [T] 



Consequently, at rotational level, the RBC experiences two distinct components of the torque. The 
first one arises from the coupling between the body motion and the fiuid vorticity, that we represent by 



the following rotational kernel 

r^(x, r,) = -7fi<5(x - R„ Q,) [a - ^(x)] = -^r5^ (^ - a;) (15) 

where 7/j is a rotational coupling coefficient and the superscript A stands for antisymmetric. This term 
depends on the body shape and orientation via the shape function and, in a linear shear flow, generates 
angular motion at constant angular velocity. 

The elongational component of the flow contributes to the orientational torque for bodies with ellipsoidal 
symmetry, being zero for spherical solutes [? ]. By deflning the stress vector f^ = cr ■ fi, where n is 
the outward normal to the surface of a macroscopic RBC, we replace the surface normal with the vector 
spanning over the entire volume of the diffused particle, fi = d6/\d6\. The associated torque is represented 
in analogy with the torque acting on macroscopic bodies [? ], by the kernel 

r^(x,ri)=a^,t"x(x-Ri) (16) 

where a is a parameter to be flxed and the superscript S is mnemonic for the symmetric contribution of 
the flow. As shown in the following, the elongational component of the torque includes an independent 
contribution arising from tank treading that will allow us to tune the parameter a based on known data 
on the tumbling to tank treading transition. 

The hydrodynamic force and torque acting on the RBC are obtained via integration over the globule 
spatial extension. Owing to the discrete nature of the mesh, the integrals are written as discrete sums, 

F, = J2 '^(^' r.) = -7t(V, - u,) (17) 

X 

Tf = E ^''(^' r^ = -M^^ - ^i) (18) 

X 
X 

where 

Ui = 5i'k\i= 2_] ^jU (20) 

X 

uii = Si-kcjj = ^5iCjj (21) 

X 

are smeared hydrodynamic flelds and the symbol -k denotes convolution over the mesh. 

The action of the forces Fj and torques Tj = Tf + Tf are counterbalanced by opposite reactions on the 
fluid side. Conservation of linear and angular momentum in the composite fluid-RBC system preserves 
the basic symmetries of the microdynamics and produces the consistent hydrodynamic response [? ] . The 
action of forces and torques on the fluid populations are expressed according to the following body force 



G = -J2[FA + ^T,xd6A 



The two exchange terms arising from the translational and rotational back-reactions produce distinct 
modifications of the fiuid velocity and vorticity. Some algebra shows that every suspended body preserves 
mass and linear momentum in the composite fiuid-RBC system since ^^Au = '^y^^p'^fpCp = — Fj. 
Similarly, every suspended body preserves the total angular momentum, since 

J2^'^ = J2Y1 ^fp^'p X (^ - R-^) = ^ 5Z (t^ X d~5?j X (x - Ri) = -Ti (22) 

X X p X 

where we have used the Lagrange rule, ax (b x c) = (a ■ c)b — (a ■ b)c and the property of the shape 



function, eq. (13). 



C. Tank Treading 

Tank Treading (TT) is the motion featured by vesicles and RBC and arises from the sliding of the 
surrounding membrane when subjected to a shearing fiow. It can be visualized by considering a material 
point anchored to the RBC membrane and moving in elliptical orbits while maintaining fixed the orienta- 
tion of the RBC with respect to the shearing direction (see Fig. fl|. The physical parameters controlling 
tank treading are the ratio of viscosity between plasma and the inner fiuid entrained within the RBC, the 
shape of the RBC and the shear rate 7 [? ]. 

In modeling blood, it is essential to include TT as it acts to orient RBC with a privileged angle with 
respect to the fiow direction. In proximity of the vessel walls, the fixed orientation induces a net lift force 
proportional to 7 that pushes the globule away from the walls [? ]. Lift forces, arising from either single- 
body or many-body effects, are thought to induce the Farhaeus-Lindqvist phenomenon, the drop of blood 
viscosity in vessels of sub-millimeter diameters, an effect with far-reaching consequences in physiology [? 

]■ 

In the general case, TT motion takes place in an apparently decoupled fashion from tumbling, as vor- 
ticity and elongational contributions have different origins. However, the motion of the rigid-body RBC is 
partially compensated by the movement of the surrounding membrane and thus effectively couples tum- 
bling and TT motion via the instantaneous orientation of the globule. At small values of the viscosity ratio, 
orientational torques prevail over the rotational ones and pure TT motion with a fixed RBC orientation 
is observed. The phase diagram as a function of the viscosity ratio provides a direct view on the RBC 
response to shear [? ]. 

An analytical treatment accounting for tank treading is due to Skalak and Keller (SK) [? ] a model 
that has enjoyed a large popularity in determining the single-body forces and torques acting on RBCs [? 
] [T]. Recently the SK model has been re-derived by Rioual et al. [? ] via a heuristic argument that we 
recall here. In 2D, when the RBC is subjected to a velocity field u = (0,7?/), the elongational component 
tends to align the RBC along the shearing principal eigenvector. For an elliptical body, the alignment has 

8 



an angle ip = j, where ip is the angle formed by the RBC largest principal axis with respect to the flow 
direction. 

In general, the torque acting on a RBC can be decomposed into three separate components. The first 
component is due to the angular motion of RBC in a quiescent fluid and produces a frictional torque given 
by —'~fRfl, where 7ij is a phenomenological coefficient that can be identified with the one introduced in 



eq.(15). The second component is due to the torque from the shearing fluid on the quiescent RBC and for 
a quiescent membrane, being related to the vorticity component of the flow. The analytical solution for the 
2D elliptical RBC without TT motion provides a total torque proportional to 7 ( ^^ ^ H — ^^-^ cos 2ip j , 
where Li^2 are two shape-dependent coefficients and where the first term of the expression is associated 
to vorticity while the second term to the elongational component of the flow. The latter provides two 
equilibrium positions for ip = ±^ and maximal torque for ip = 0,^. 

Finally, the effect of tank treading is determined by considering the local frame moving together with 
the material point at velocity V^t and with the infinitesimal membrane element experiencing a force 
dFj"T oc Vtt and torque Ttt = § (IYtt x (x — R) . Consequently, TT couples to both the rotational 
and elongational flow components and results in a net torque Ttt oc cos(2'?/'), that is, a torque enjoying 
the same angular symmetry of the mechanism associated to the rigid body response. Given the three- 
dimensional, diffused nature of the numerical model, the actual presence of the cellular membrane is not 
explicitly considered. Instead, the effect of TT is controlled by tuning the intensity of the elongational 



torque via the adimensional factor a of eq. (16). 



D. Viscosity contrast 

RBC are carriers of an internal fluid composed of several hundred millions of hemoglobin proteins and 
being much more viscous than plasma. The inner fluid contributes significantly to the dissipation of energy 
as mediated by the shearing modes activated by the RBC membrane. The consequence is the steep raise 
in the apparent viscosity with the hematocrit level. We incorporate the effect of the viscosity contrast by 
considering a local enhancement of the LB fluid viscosity within the RBC shape according to the following 
BGK relaxation time 

T{x) = % + ^hY,Oi (23) 

i 

where To corresponds to the viscosity of pure plasma, A is viscosity enhancement factor, and 9i is a smooth 
version of the ellipsoidal characteristic function. The latter is represented as ^j = 1 — (1 — ^j)'' where the 
parameter k governs the smooth transition between the inner and outer fluids, chosen as /t = 20. By 
choosing z/q = 1/6 and A = 2, the ratio between inner {6i ~ 1) and outer {6.1 ~ 0) viscosities is equal to 5. 



E. Excluded Volume Interactions 

In order to account for the RBC-RBC repulsive forces, soft-core mechanical forces and torques are given 
by the Gay-Berne (GB) potential [? ], as widely used in the liquid crystals community. The GB potential 
inhibits interpenetration of pairs of RBCs by introducing an orientation-dependent repulsive interaction 
derived from the Lennard- Jones potential {(pLjiRij) = ^[(^/-Rji)^^ ~ i'^/Rij)^]^ where e is the energy scale 
and a the "contact" length scale). The GB potential extends the spherically symmetric Lennard- Jones 
potential to ellipsoidal symmetry, where the potential depends if two RBCs have mutual orientation as 
face-to-face (maximal repulsion), side-to-side (minimal repulsion) or an arbitrary orientation between the 
two bodies (intermediate case). 

In the following, we shall employ the form of the GB potential that handles the interaction between 
particles of different eccentricity, such as a mixture of ellipsoidal and spherical particles. This flexibility 
allows to handle biofluids composed of particles with different shapes by employing the same analytical 
form of the potential. 

Given the principal axes (01,1,01,2)^1,3) of the i-th globule, the ellipsoidal shape associated to the ex- 
cluded volume interactions is constructed according to the shape matrix Si = diag(ai,i, 0^,2, 01,3) and the 
transformed matrix Aj = QjSfQf in the laboratory frame. The pair of particles i,j at distance Rjj 
experiences a characteristic exclusion distance Cij that depends on the globule-globule distance, shape and 
mutual orientation, written as 

^., = -^ (24) 

where the matrix H^- = Aj + Aj has been introduced. 

A purely repulsive exclusion potential is given by the pairwise form [? ? ] 

1 4>2 

r ''J ^mm ^ 1 

ij 

where eo is the energy scale and o"^*" is a constant, both parameters being independent on the ellipsoidal 
mutual orientation and distance. For two identical oblate ellipsoids, o"^*" corresponds to a contact distance 
of the two particles having face-to-face orientation. In general, by considering the minimum particle 
dimension a™" = min(aj,i, aj,2, Oj.s) then 



with 






(a™")^ + (a™") (2^ 

10 



The mutual forces and torques exerted between the i,j pair are computed accordingly, as described in the 
appendix for the sake of completeness. 

We note here two aspects. First, the softness of the RBC can be improved by taking small values 
for the energy scale eo. Clearly, this type of softness does not arise from hydrodynamic forces but from 
RBC-RBC mechanical forces, and emulates the good packing attitude of RBCs at high hematocrit levels. 
Secondly, the ellipsoidal exclusion shape can be modified to reproduce more closely the biconcave discoid 
via the Bates-Lockhurst version of the GB potential [? ]. At present, we take the ellipsoidal shape to be 
representative of RBC, while further refinements can be trivially included. 

When simulating blood fiows in a confining environment, we associate to each wall mesh point a spherical 
particle that acts as a repulsive center for the ellipsoidal RBC and prevents the globule from leaking out 
of the confining vessel, under the action of the same GB potential. The wall-RBC characteristic distance 
'^wRBC i^ chosen to be half the mesh spacing ^. 

F. Summary of the computational model 

To recap, the basic assumption of the model is to handle a suspension made of plasma and RBC. 
The former is a Newtonian fiuid that reproduces Navier-Stokes dynamics at virtually arbitrary Reynolds 
numbers and in arbitrary geometries. The physiological conditions of Reynolds < 2000 and shear rates 
< 500 s~^ can be accessed in simulation without posing limitations in terms of numerical robustness and 
feasibility. The model for RBCs has the following characteristics: 

• it is based on the body forces and torques acting on the RBC rather than on surface forces by a 
proper decomposition of the translational, tumbling and orientational components of the fiow. A 
similar scenario is found in Faxen's law for zero- Reynolds fiows [? ]. Our strategy does not track 
explicitly the RBC membrane or handle in anyway the internal cytoskeleton; 

• RBCs are active scalars with hydrodynamic shape that is fixed and ellipsoidal. The eccentricity of 
the RBC can be tuned, however in the following we choose to work with C,x = ^ and ^y = C,z = 2, 

47r{S'/47r)a/2 



corresponding to volume V ^ 134, surface S ^ 139 and reduced volume v = ,^^X .3,2 ~ 0.87, that 



should be compared with v = 0.65 for human RBCs; 

The near-field hydrodynamic interactions are well captured by the model, in particular as regarding to 
translational and orientational mobilities and long-time tails arising from hydrodynamic interactions. 
Shape fiuctuations can be easily introduced, e.g. by means of the Maffettone-Minale extension [? ], 
by allowing volume-preserving deformations while still maintaining the ellipsoidal symmetry; 

the viscosity contrast is introduced by allowing for a modification of the relaxation time within 
the RBC region. In practice, this means solving the same LB algorithm throughout the system, 

11 



but with a local relaxation time and thus, higher dissipation within the RBC region. In terms of 
computational feasibility, a wide range of viscosity contrast can be explored, but in the following we 
will focus on a given viscosity contrast as pertaining to RBCs; 

• interpenetration between RBCs is avoided by Gay-Berne forces and torques. The degree of softness 
of particles can be tuned, as much as some attracting component of the interaction due to adhesive 
forces or the presence of fibrinogen, can be easily accommodated within the model. 

In summary, the main assumption of the model, the neglect of shape fluctuations, is largely compensated 
by its main strength, the possibility to handle large-scale systems of physiological relevance with state-of- 
the-art computer hardware. This is a major point of the model that will not be addressed in detail here, 
but has been discussed in refs.[? ? ? ]. 

III. NUMERICAL RESULTS 

The frictional response of a single RBC in plasma is first analyzed by computing the Stokes response 
of an oblate ellipsoidal particle having two possible orientations with respect to the motion direction. For 
translational displacement, these are the frontal and side-wise motions. For the rotational motion, these 
correspond to spinning around two principal directions corresponding to the smallest and largest radii. 
For this preliminary test, we deactivate the modulation of viscosity by the presence of the RBC, that 
is, we take the same viscosity for plasma and the entrained fluid without distinction and throughout the 
simulation domain. 

As shown in ref. [? ] for a model of suspended point-like particles, the effective mobility is given by 
the sum of two components, the mobility associated to the bare frictional parameters, 75- and 7/j, and the 
effect of the hydrodynamic field induced on the surrounding solvent that sustains the motion by increasing 
the particle roto-translational mobilities. In addition, the hydrodynamic components to mobility contains 
a Stokes-like component that is renormalized by the presence of the numerical finite-spacing mesh. We 
expect to recover the same qualitative behavior for the diffused ellipsoidal particles. 

Fig|2] shows the computed particle mobilities as a function of the frictional parameters 7^? and 7^ 
for a RBC of mass M = 10 and inertia Ix,y,z = 1000 (all data are expressed in lattice units unless 
otherwise stated). The particle is embedded in a simulation box of size 60 x 60 x 60 containing a quiescent 
fluid. At infinite friction, the intercepts correspond to the mesh-induced spurious frictional forces. For 
translational motion, the intercepts correspond to two mobilities, /i = l^^^^, = 6.6 and 10.3, for frontal 
and side-wise translation motion, respectively. By associating a residual hydrodynamic radius a^'^^^ = 
—^ — = 0.06 and 0.03 to frontal and side-wise motion, respectively, it is apparent that the residual 
radius is directly proportional to the size factor ^q, governing the shape function. From these data it is 
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Figure 2: Translational and rotational mobilities as a function of the coupling parameters 7^ and 7/?. Circles 
correspond to frontal (filled symbols) and lateral translation (open symbols). Squares are for rotations around the 
smallest principal radius (filled symbols) and around the largest principal radius (open symbols). The lines are 
guides for the eye. 

clear that the residual Stokes radii are much smaller than the mesh spacing (being Ax = 1 in lattice units) 
and modulating the hydrodynamic response to achieve a bulkier suspended body requires increasing the 
coupling parameter to 77- > 10. High values of 7^ spoil the numerical robustness of the method. In fact, 
a simple estimate shows that numerical stability of the LB method is for body forces \F\ < 0.1, that is 
It ^ 0.1. The translational data show that the coupling method alone reproduces the far-field Stokes 
response but is unable to modulate the viscosity to significant values with the hematocrit level, unless a 



different mechanism is introduced to improve the internal dissipation, such as the one encoded by eq. (23). 
For angular motion, the situation is slightly different since we find the intercepts 7™'='''^ = ^ At~^ 
and ^ At^^ for rotations around the smallest and largest principal radii, respectively. The intercepts 



38 



correspond to the residual rotational radii a^ 



mesh 



h\ 1/3 



Sirup 



1.84 Ax and 1.75 Ax, respectively. In the 



angular motion, the mesh-induced friction is larger than the translational counterpart (in fact, comparable 
to the mesh spacing) and exhibits a weak dependence on the direction of spinning direction, so that it can 
be considered independent on the latter. 



The effect of tank treading is analyzed by placing a RBC in a shear field generated according to the 
method of ref. [? ]. Fig. IS] shows the time dependence of the elongational component of the torque for 
a single RBC in a shear flow. Different values of the enhancement factor a are reported having chosen 
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Figure 3: Symmetric component of the torque of a single RBC in a shear field and for different values of the tank 
treading parameter a. 

It = iR = 10^^. It is found that the value a = 0.063 corresponds to the transition to pure TT motion, 
where the elongational and vorticity-related torques exactly compensate, generating a constant orientation 
with angle ~ 40° as a fixed point of the dynamics. For a bounded fluid in a two-dimensional slab the flxed 
angle is found to be in the range 20° — 35°, depending on the viscosity contrast and flow curvature, with 
a lift force being linearly proportional to the local shear rate and to y^"^, where y is the distance of the 
RBC from the wall, in agreement with experimental and theoretical results [? ? ]. 

Overall, the torques generate a well-behaved rotational motion and in line with the expected trend. In 
the following we arbitrarily set the tunable parameter a = 0.04, that is, slightly below the transition to 
tank treading, and proceed further with tests regarding the rheological response of a dense suspension. 
The tests in the bounded fluid conflrm the tumbling regime for the chosen RBC shape and viscosity ratio. 
It should be borne in mind that the balance between tumbling and tank treading in a dilute suspension 
critically affects the lift forces and thus the viscosity of the suspension [? ]. In this respect, the flxed-shape 
approximation of the present model produces hydrodynamic forces that are qualitatively correct. In the 
future, the model can be tuned more flnely by optimizing the shape and viscosity contrast. 



The viscosity of the suspension is computed by using a cylindrical channel of radius 50 fim and for 
different hematocrit levels. In Fig. |4] we report the relative viscosity {i^rei = Vapp/Vo) where i^app is the 
apparent viscosity measured in the channel at flnite hematocrit and rjo the viscosity in the same channel 
at zero hematocrit. Since t] oc Q, where Q is the volumetric flow rate, rirei is computed as the ratio of flow 
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Figure 4: Relative viscosity in a channel of radius 50 ^m for different hematocrit levels as compared to the 
experimental data of Pries et al. [? ] (solid curve). Data are for the enhanced dissipation mechanism of eq. (23) 
with A = 2 (circles) and without enhancement (A = 0) (squares). 

rates at finite and zero hematocrit. Fig. |4]also reports the data on viscosity by setting the enhancement 



factor of eq. (23) to A = 0. The latter produce a weak modulation of viscosity with hematocrit, while the 



data with A = 2 exhibit an excellent agreement with the experimental results of ref. [? ]. 



A crucial phenomenology of blood circulation is the decrease of viscosity as the vessel radius falls below 
100 fim, namely, the Farhaeus-Lindqvist effect [? ]. This effect is usually ascribed to the formation of a 
cell-free layer in proximity of the vessel walls. The origin of such depletion is still uncertain but the lateral 
forces that push the RBCs away from the vessel walls are retained to have different causes, such as tank 
treading and cell deformation [? ] , adhesive properties of RBC or shear-induced migration. In the current 
version of our model, we do not probe the effects of cell deformation. However, the simulation reveals a 
distinct RBC depletion in proximity of the walls, as shown in Fig. |5| We report the variation of the cell-free 
layer with the hematocrit and vessel diameter as compared to the experimental data of Bugliarello and 
Sevilla [? ]. The numerical results reproduce the experimental data quite well, lending good confidence 
in the numerical model at vessel diameters below the 100 fim radius. The typical cell-free layer observed 
in simulation for two channels of different radii is illustrated in Fig. [5} inset. The progressive depletion 
of RBC as compared to the plasma content is the leading cause of fluidization of the suspension at small 
channel radius. 
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Figure 5: Size dependence of the cell-free layer with the vessel diameter and hematocrit level of 10% (diamonds), 
20% (squares) and 50% (circles), as compared to the experimental data of Bugliarello and Sevilla [? ] (star 
symbols). The lines are guides for the eye. Inset: Radial density profiles of RBC, for R = 10 fim (solid line) and 
R = 20 nm (dashed line), illustrating the cell-free layers in proximity of the vessel wall. 

The non-Newtonian behavior of the suspension is further exhibited by the velocity profiles of RBC for 
different hematocrit levels and vessel diameters, as shown in Fig. |6| As the hematocrit level increases, the 
Poiseuille-like parabola modifies into flatter profiles next to the vessel centerline and in a large extension 
of the channel, whereas in proximity of the walls, the profiles have large slopes and strong dissipation, in 
particular for the narrower vessels. We further notice a good match of the plasma and RBC velocities for 
all vessel radii, a matching that is typically lost as the vessel radius becomes comparable to the RBC size 
(data not shown). 



In applying the model to physiological conditions, we consider a realistic bifurcating vessel at 50% 
hematocrit level, as depicted in Fig. [7j The bifurcation is extracted as part of a coronary arterial system 
[? ], and is made of a parent vessel of radius ~ 100 /im and one daughter branch having approximately 
the same size of the parent vessel, and a second daughter branch of radius ~ 80 fim. 

The flow is induced by imposing parabolic flow conditions at the inlet and constant pressure at the 
outlets, with the method described in [? ], producing a shear rate of 80 s^^ averaged over the entire 
bifurcation at steady state, and reaching as high as 150 s^^ in proximity of the arterial wall, values in line 
with the typical shear rates observed in human coronary arteries, typically in the 30 — 450 s~^ range [? 
]. As the snapshot reveals, RBCs organize in several different ways throughout the bifurcation and also 
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Figure 6: Velocity profiles for vessel radius of 10 (upper panel), 25 (mid panel) and 50 fim (lower panel). Data 
correspond to hematocrit levels of 35% (solid lines) and 0% (dashed lines). 

depending on the value of the shear rate (data not shown). The local organization of RBCs in rouleaux 
is clearly visible, the typical stack often observed in static or flow conditions, and mostly destroyed as the 
shear rate increases. 



The flow in the bifurcating channel for different values of the shear rate is analyzed and the shear 
thinning behavior is illustrated in Fig. [8| where viscosity decreases with the shear rate until the Newtonian 
plateau is reached. 

As expected, shear thinning correlates with the formation of rouleaux, as seen by inspecting the sta- 
tistical distribution of cluster size. The histogram of the rouleaux size is reported in Fig. [8} inset. As 
shear rate increases, fewer and smaller structures are detected. However, a central observation is that the 
local structuring of RBCs strongly depends on the morphological details of the vessel. In fact, we notice 
that in proximity of the bifurcation a shoulder in a daughter vessel creates smaller flow velocity and some 
stagnation of RBCs. In correspondence of this region, rouleaux show resilient character even for shear 
rates above 50 s^"*^. 

The uneven distribution of RBCs, with the attendant stagnation and persistence of rouleaux in speciflc 
regions, can have signiflcant impact on the distribution of shear stress. The latter quantity strongly cor- 
relates with the formation of plaques in the endothelium (atherogenesis) and the subsequent development 
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Figure 7: Snapshot of the bifurcating vessel showing the organization of RBCs at hematocrit of 50 % and average 
shear rate of 80 s~^. 




shear rate (s ) 



Figure 8: Shear thinning of the suspension in the bifurcating channel for different values of the shear rate. Inset: 
Distribution of rouleaux as a function the cluster size S, for shear rate of 8 s^^ (solid bars) and 16 s^^ (open bars). 

of atherosclerosis [? ]. In particular, low levels of shear stress, as due to disturbed flow patterns and 
stagnation regions, trigger the growth of plaques. It is therefore of great relevance to inspect the role of 
RBCs in modulating the shear stress distribution on the wall of the vessel under study. 

For this analysis, we employ the method described in ref. [? ] and compute the endothelial shear 
stress as time averages over the evolution of the plasma-RBC composite system. Fig. [9] illustrates the 
distribution of shear stress for different hematocrit levels. The plot reveals the strong effect of the RBCs 
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Figure 9: Time-averages of shear stress in the bifurcating channel for pure plasma (left panel), hematocrit level of 
35% (mid panel) and 50% (right panel). The flow corresponds to an average shear rate of 80 s^^. 

throughout the system and in particular the great fluctuations in one daughter vesseL While the overall 
shear stress distribution is somehow preserved at different hematocrit levels, important local modiflcations 
are induced by RBCs. In particular, in proximity of the vessel shoulder, the RBC structuring induces 
smaller values of the shear stress, followed by larger values next to the inner side of the bifurcation. 
Given these observations, it is expected that the particulate nature of blood has signiflcant effects in 
the identiflcation and prediction of vulnerable regions of the large-size vascular system by computational 
methods. 



IV. CONCLUSIONS 



In this paper, a new model for suspended red blood cells has been introduced with the goal of studying 
large-scale hemodynamics in cardiovascular systems. The underlying idea is to represent the different 
responses of the suspended bodies, emerging from the rigid-body as much as the vesicular nature of the 
globule, by distinct coupling mechanisms. These mechanisms are entirely handled at kinetic level, that 
is, the dynamics of plasma and red blood cells is governed by appropriate collisional terms that avoid 
to compute hydrodynamic forces and torques via the Green's function method, as employed in Stokesian 
dynamics [? ] . The fundamental advantage of hydrokinetic modeling is to avoid such expensive route and, 
at the same, enabling to handle flnite Reynolds conditions and complex boundaries or irregular vessels 
within the simple collisional approach. At macroscopic scale, the non-trivial rheological response emerges 
spontaneously as a result of the underlying microdynamics. 
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The second, crucial advantage of the proposed method is to enable the investigation of systems of 
physiological relevance whose phenomenology typically spans several decades in time and space. To this 
purpose, massive parallel computers must be deployed by exploiting a numerical method that scales over 
thousands of processors and above. To a large extent, both the Lattice Boltzmann method and the red 
blood cell dynamics have been proved to scale over traditional CPU-based computers such as on Blue 
Gene architectures, as much as over massive assemblies of Graphic Processing Units (see ref. [? ] and 
references therein). 

The numerical results presented in this paper have shown that the particulate nature of blood cannot be 
omitted when studying the non-trivial rheology of the biofluid and the shear stress distribution in complex 
geometries. Regions of low shear stress can appear as the hematocrit reaches physiological levels as a result 
of the non-trivial organization of red blood cells and the irregular morphology of vessels. These results 
cannot can have far reaching consequences in real-life cardiovascular applications, where the organization 
of red blood cells impacts both the local flow patterns and the large-scale flow distribution in complex 
vascular networks. 
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Appendix A: Rigid Body Dynamics 

We report here the numerical solution of the roto-translational dynamics for RBCs. The rotational 
dynamics is customarily solved in the body frame where the equation relating angular momentum tt and 
angular velocity f2 is 

n = in (Al) 

where, for ease of notation, we have omitted the RBC index. The rigid-body dynamics is given by 

f2 = I-i [(m) X fl] + I- V (A2) 

and the rotational matrix obeys the equation 

Q = Qskew{n) (A3) 

where, for a vector with components v = {vn,Vt,Vg), skew{\) denotes 

^ Vg -Vt^ 

skewiy) = -Vg Vn (A4) 

(^ Vt -Vn J 
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In absence of hydrodynamic forces and torques, the dynamics is Hamiltonian and the updating algo- 
rithm is given by the scheme described in ref. [? ]. For 7t = 7i? = the algorithm is symplectic and 
provides excellent conservation of energy. In presence of hydrodynamic forces and torques, the algorithm 
is generalized as in the following steps: 

1. Propagate linear and angular momenta for half-step according to the combination of mechanical and 
hydrodynamic forces F'^"™^ and torques r^""^^ , 

nh/2 = no + h-'r^r' (A5) 



h 

2M' 



V,/2 = Vo + ;t77Fo°™' (A6) 



where 



^c^omb ^ F™'^"'' + 7rMiio (A7) 

tT' = tZ'^ + ^rIG^o (A8) 

2. Update linear and angular momenta for half-step according to the frictional forces 

^h/2 = e-'^^/^no (A9) 

3. Propagate positions as 

Rh = Ro + Wv2 (All) 

4. Update the angular state by taking the rotation matrix TZ{a,P,'~f), associated to rotation around 
the axes n, t and g with angles a, (5 and 7, respectively. Propagate the rotational state according 
to the sequence of five sub-rotations 



(a) 7^ = 7^(^,0,0) 



(b) 7^ = 7^(o,f,o) 



(c) 7^ = 7^(o, 0, /i7) 



(d) 7^ = 7^(o,f,o) 




Q = Q7^t 
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(e) 7^ = 7^(f,0,0) 




5. Compute the mechanical forces and torques from the excluded volume interactions, F^ and r^ 



.mech 



and the fluid velocity and vorticity via Eqs.(20 21), and construct the combinations 



_^„,,„„ ^mecn , „, t,~ 



h — ^' h 

_comb _Tnech 



(A12) 
(A13) 



6. Update linear and angular momenta for half-step according to the frictional forces 



(A14) 
(A15) 



7. Finalize the linear and angular momenta for half-step according to combined forces and torques 

h 






VV2 + 



2M 



^comb 
h 



h. 



l_,comb 






(A16) 
(A17) 



Appendix B: Gay-Berne forces and torques 



The calculation of mechanical forces and torques according to the Gay-Berne potential, eq. (26), is 
reported here. For the interacting pair i,j , the pairwise mechanical force is given by 



^GB 



du 



«j 



dRij 



1 du 



I] '^Yi] 



Rijd(t),jd%j ^^ '' '' 



duij_^ 



1 duij 



'^ij x'^ij ' ^^ijj^^ij 



where the following vector has been defined 



l^ij — tlij ■ t^ij 



The terms appearing in eq. (Bl) are computed as 

duj. 



OR, 



-24en 



■ij 



(2pr-p^) 



and 



dui 



V 



-24e, 






2a,7" 



(Bl) 



(B2) 



(B3) 



(B4) 



The torques are calculated according to the prescription in ref. [? ] and transformed back to the body 
frame 

rf^ = QnRc.xf.J (B5) 
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and 



where 



r™ = Qj(R,,xf,._,) (B6) 



' *«-=« (B7) 



Rij d(f)ij 
and 

It is easy to show that torques obey conservation of angular nionientuni since rf^ + r^^ = — Rjj x F^^. 



[1] In its original formulation, the Skalak-Keller model is two-dimensional and considers the torques acting on 
elliptical RBC in Stokes flow and obeying the superposition principle, that is, with tumbling and tank treading 
taken as independent motions. Moreover, the RBC is considered as passive scalar with no back-reaction on the 
fluid. The SK model predicts tumbling to TT transition independent on the shear rate 7. 
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